Introduction à la manipulation de données géospatiales avec R

INTRODUCTION

The two most important R packages for manipulating spatial objects are:

sf
vectors, points, lines, polygons, shapefiles and geopackages

terra
rasters or regularly gridded data


You may see some code or package that still rely on the packages sp, rgdal, rgeos or raster, but these are now mostly retired and were entirely replaced by sf and terra.

Another important package is stars for spatiotemporal arrays and raster and vector data cubes, but we won’t talk about this package in this workshop.

VECTOR DATA (sf)

Reading

First, let’s look in the folder where the shapefile is located.

list.files("/home/frousseu/Documents/github/spatialR/data")

##  [1] "carreteras.dbf"        "carreteras.prj"        "carreteras.sbn"       
##  [4] "carreteras.sbx"        "carreteras.shp"        "carreteras.shp.xml"   
##  [7] "carreteras.shx"        "cat.txt"               "cdem_dem_021E_tif"    
## [10] "cdem_dem_021E.tif"     "cdem_dem_031F_tif"     "cdem_dem_031F_tif.zip"
## [13] "cdem_dem_031F.tif"     "cdem_dem_031G_tif"     "cdem_dem_031G_tif.zip"
## [16] "cdem_dem_031G.tif"     "cdem_dem_031H_tif"     "cdem_dem_031H.tif"    
## [19] "rast.tif"              "roads.dbf"             "roads.gpkg"           
## [22] "roads.prj"             "roads.shp"             "roads.shx"            
## [25] "test.dbf"              "test.gpkg"             "test.prj"             
## [28] "test.shp"              "test.shx"

Reading shapefiles is done with the function st_read.

library(sf)

roads <- st_read("/home/frousseu/Documents/github/spatialR/data/carreteras.shp")

## Reading layer `carreteras' from data source 
##   `/home/frousseu/Documents/github/spatialR/data/carreteras.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1171 features and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 123569.7 ymin: 1963648 xmax: 309483.2 ymax: 2128266
## Projected CRS: WGS 84 / UTM zone 16N

Plot it

This object can be plotted simply with the function plot which has a method for sf objects. The st_geometry is used to only plot the geometry of the object, but not the attributes.

plot(st_geometry(roads), axes = TRUE)

What exactly is the object roads?

class(roads)

## [1] "sf"         "data.frame"

head(roads)

## Simple feature collection with 6 features and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 269922.2 ymin: 1964531 xmax: 292506.3 ymax: 1976234
## Projected CRS: WGS 84 / UTM zone 16N
##   CLASIF   LENGTH DESCRIPCIO      PAíS                       geometry
## 1     C4  747.849    Veredas Guatemala LINESTRING (270613.6 196999...
## 2     C4 1303.025    Veredas    Belice LINESTRING (281424.3 197242...
## 3     C4 2055.814    Veredas    Belice LINESTRING (292506.3 196456...
## 4     C4 1807.620    Veredas    Belice LINESTRING (277670.8 196523...
## 5     C4 2661.517    Veredas    Belice LINESTRING (272900.7 197484...
## 6     C4  547.114    Veredas    Mexico LINESTRING (270909.1 197623...

It is mostly a data.frame to which we add a geometry column.

Almost just a data.frame

head(as.data.frame(roads))

##   CLASIF   LENGTH DESCRIPCIO      PAíS                       geometry
## 1     C4  747.849    Veredas Guatemala LINESTRING (270613.6 196999...
## 2     C4 1303.025    Veredas    Belice LINESTRING (281424.3 197242...
## 3     C4 2055.814    Veredas    Belice LINESTRING (292506.3 196456...
## 4     C4 1807.620    Veredas    Belice LINESTRING (277670.8 196523...
## 5     C4 2661.517    Veredas    Belice LINESTRING (272900.7 197484...
## 6     C4  547.114    Veredas    Mexico LINESTRING (270909.1 197623...

Plot a single attribute

plot(roads["CLASIF"], lwd = 4, axes = TRUE)

Plot (all) attributes

plot(roads, lwd = 4, axes = TRUE)

Plotting only the geometry

Often, we only want to plot the geometry to see if everything looks OK. This is done with the st_geometry function.

plot(st_geometry(roads), lwd = 4, axes = TRUE)

Interactive map

Alternatively, you can use the mapview package for a more interactive plot.

library(mapview)
mapview(roads)

Writing

We can write in different formats. The geopackage format (.gpkg) is used more and more compared to the .shp format.

# As a shapefile
st_write(roads, dsn = "/home/frousseu/Documents/github/spatialR/data/roads.shp", append = FALSE)  # append = FALSE overwrites

## Deleting layer `roads' using driver `ESRI Shapefile'
## Writing layer `roads' to data source 
##   `/home/frousseu/Documents/github/spatialR/data/roads.shp' using driver `ESRI Shapefile'
## Writing 1171 features with 4 fields and geometry type Line String.

# As a geopackage
st_write(roads, dsn = "/home/frousseu/Documents/github/spatialR/data/roads.gpkg", append = FALSE)

## Deleting layer `roads' using driver `GPKG'
## Writing layer `roads' to data source 
##   `/home/frousseu/Documents/github/spatialR/data/roads.gpkg' using driver `GPKG'
## Writing 1171 features with 4 fields and geometry type Line String.

list.files("/home/frousseu/Documents/github/spatialR/data")

##  [1] "carreteras.dbf"        "carreteras.prj"        "carreteras.sbn"       
##  [4] "carreteras.sbx"        "carreteras.shp"        "carreteras.shp.xml"   
##  [7] "carreteras.shx"        "cat.txt"               "cdem_dem_021E_tif"    
## [10] "cdem_dem_021E.tif"     "cdem_dem_031F_tif"     "cdem_dem_031F_tif.zip"
## [13] "cdem_dem_031F.tif"     "cdem_dem_031G_tif"     "cdem_dem_031G_tif.zip"
## [16] "cdem_dem_031G.tif"     "cdem_dem_031H_tif"     "cdem_dem_031H.tif"    
## [19] "rast.tif"              "roads.dbf"             "roads.gpkg"           
## [22] "roads.prj"             "roads.shp"             "roads.shx"            
## [25] "test.dbf"              "test.gpkg"             "test.prj"             
## [28] "test.shp"              "test.shx"

Building shapefiles

cat <- read.table("/home/frousseu/Documents/github/spatialR/data/cat.txt", header = TRUE)
head(cat)

##           X        Y Attack
## 1 -89.34188 18.49952      0
## 2 -89.41309 18.47991      0
## 3 -89.26248 18.60019      1
## 4 -89.19920 18.59252      1
## 5 -89.30334 18.59921      1
## 6 -89.23133 18.59242      1

This is a simple data.frame with X and Y columns representing longitudes and latitudes. The column Attack contains whether a location was concerned by a jaguar attack or not.

To transform this data.frame to a spatial object, we just have to do:

cat <- st_as_sf(cat, coords = c("X", "Y"), crs = 4326)

The object returned is now an sf data.frame. The code 4326 is a special id for the lat/lon WGS84 coordinate reference system (crs).

Plot it

Here is what it looks like.

plot(st_geometry(cat), col = adjustcolor(ifelse(cat$Attack == 0, "blue", "red"), 0.4), pch = 16,
    cex = 1.5, axes = TRUE)

Add the 2 layers in an interactive map

mapview(roads) + mapview(cat)

– Exercice 1 –

Build a spatial data.frame with 100 random locations centered on the CEF conference, plot the results with R and mapview and write the results to a shapefile or a geopackage (.gpkg). Read the shapefile you wrote to make sure everything worked well.

Hint: use the function runif or rnorm to generate random values within a given range. For example:

lat <- rnorm(100, mean = 45, sd = 0.1)  # 100 random values with a mean of 45 and an sd of 0.1
lon <- rnorm(100, mean = 72, sd = 0.1)

COORDINATE REFERENCE SYSTEM (CRS)

quebec <- st_read("/home/frousseu/Documents/github/spatialR/québec.gpkg")

## Reading layer `québec' from data source 
##   `/home/frousseu/Documents/github/spatialR/québec.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -80.27292 ymin: 44.9912 xmax: -57.10779 ymax: 63.68542
## Geodetic CRS:  WGS 84

st_crs(quebec)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

Different CRS

par(mfrow = c(1, 2))
plot(st_geometry(quebec), axes = TRUE)
quebec <- st_transform(quebec, 6622)
plot(st_geometry(quebec), axes = TRUE)

As always, when playing with spatial data, one needs to be aware of coordinate reference systems (CRS). Assigning and transforming CRS in sf can be done with the st_transform function.

Red Trillium

First, let’s get some data from GBIF using the rgbif package. We will get occurrence data for a spring plant species, the Red Trillium (Trille rouge, Trillium erectum).

Get occurrences from GBIF

library(rgbif)

occs <- occ_search(scientificName = "Trillium erectum", hasCoordinate = TRUE, stateProvince = "Québec",
    country = "CA", limit = 5000)
trille <- st_as_sf(occs$data, coords = c("decimalLongitude", "decimalLatitude"))
plot(st_geometry(trille), axes = TRUE)

st_crs(trille)

## Coordinate Reference System: NA

We haven’t assigned a crs to the sf object.

EPSG and CRS

See https://epsg.io/ to find CRS EPSG codes.

4326 - WGS 84 – WGS84 - World Geodetic System 1984, used in GPS
32198 - NAD83 / Quebec Lambert
26918 - NAD83 / UTM zone 18N
9820 - Lambert Azimuthal Equal Area
6622 - NAD83(CSRS)v2 / Quebec Lambert

etc.

Assign a CRS

Assigning a coordinate reference system can be done with the functions st_crs or it can be done when transforming a data.frame to a spatial data.frame with st_as_sf.

st_crs(trille) <- 4326
st_crs(trille)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

With st_as_sf

trille <- st_as_sf(occs$data, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)
plot(st_geometry(trille), axes = TRUE)

Change a CRS

# 32198 NAD83 / Quebec Lambert
trille <- st_transform(trille, 32198)
plot(st_geometry(trille), axes = TRUE)

Use a common CRS

When plotting or working with multiple spatial objects, the same CRS needs to be used!

plot(st_geometry(quebec))
plot(st_geometry(trille), add = TRUE)

A note on geodata and rmapshaper

The geodata package returns a SpatVector which we transform to a sf object.

The ms_simplify function from the rmapshaper package is very useful to reduce the precision of polygons and maintaining the topology (and adjacent polygons).

library(geodata)

canada <- gadm("GADM", country = "CAN", level = 1) |>
    st_as_sf()  # transform a SpatVector to an sf data.frame


library(rmapshaper)

canada <- ms_simplify(canada, 0.001)  # simplify shape and keep 0.1% of points

Multiple CRS

Let’s use Canada to more easily see the differences between some CRS.

crs <- c(4326, 2958, 32198, 3995)

par(mfrow = c(2, 2), mar = c(2, 3, 2, 0))

for (i in 1:length(crs)) {
    canada2 <- st_transform(canada, crs = crs[i])
    plot(st_geometry(canada2), axes = TRUE, main = paste("epsg", crs[i]), cex.main = 0.7)
}

Geographic vs. Projected CRS

Geographic coordinates are coordinates on a sphere or on an ellipse, while projected coordinates are defined on a flat 2D surface. Usually, geographic coordinates are in latitudes/longitudes and projected coordinates are in meters. As seen earlier, this is important when doing certain spatial operations. All functions from package rgeos require that spatial objects are projected. Otherwise, functions often do not work properly because they assume that calculations are done on cartesian coordinates.

st_bbox(canada)

##       xmin       ymin       xmax       ymax 
## -141.00687   41.67693  -52.63496   83.10833

st_bbox(st_transform(canada, crs[3]))

##       xmin       ymin       xmax       ymax 
## -3754088.6  -146642.3  1182482.0  4557090.5

CRS and EPSG

The EPSG code 4326 is equivalent to giving the "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" CRS to the data.

– Exercice 2 –

Get some occurrences for the monarch butterfly (Danaus plexippus) using rgbif and the occ_search function and plot them for Quebec using the gadm function from the geodata package. Use the UTM zone 18N CRS.

Spatial operations

Most spatial operations are available in the package sf. One of the requirement for functions to work well for spatial operations in sf is that objects are in a projected coordinate system (i.e. not in lat/lon 4326). This can be checked with st_is_longlat.

st_is_longlat(trille)

## [1] FALSE

st_is_longlat(quebec)

## [1] FALSE

st_is_longlat(canada)

## [1] TRUE

Intersections

The function over from package sp is used to determine whether different entities overlap. It returns a vector or a data.frame with the identities or the characteristics of overlapping elements.

outaouais <- st_read("/home/frousseu/Documents/github/spatialR/outaouais.gpkg")

## Reading layer `outaouais' from data source 
##   `/home/frousseu/Documents/github/spatialR/outaouais.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -75.90876 ymin: 45.3727 xmax: -75.34527 ymax: 45.59912
## Geodetic CRS:  WGS 84

trille <- st_transform(trille, st_crs(outaouais))

plot(st_geometry(outaouais), axes = TRUE)
plot(st_geometry(trille), add = TRUE)

Interactive map

mapview(outaouais) + mapview(trille)

Intersecting

s <- st_intersects(trille, outaouais) |>
    lengths()

trilleout <- trille[s == 0, ]
trillein <- trille[s > 0, ]

plot(st_geometry(outaouais), axes = TRUE)
plot(st_geometry(trillein), pch = 16, col = "red", add = TRUE)
plot(st_geometry(trilleout), pch = 16, col = "blue", add = TRUE)

Intersecting 2nd version

Another more compact way of getting only the occurrences that are in the outaouais polygon.

trillein <- trille[outaouais, ]

Operations on points

Here are some common operations that can be done on points (but also on lines or on polygons).

x <- trillein
st_is_longlat(x)

## [1] TRUE

x <- st_transform(x, 32198)  # use a projected crs

par(mfrow = c(2, 3), mar = c(0, 0, 1, 0))

x2 <- st_union(x) |>
    st_centroid() |>
    st_geometry()

plot(st_geometry(x), main = "st_centroid")
plot(x2, pch = 16, cex = 2, col = "red", add = TRUE)

x2 <- st_union(x) |>
    st_convex_hull() |>
    st_geometry()

plot(st_geometry(x), main = "st_convex_hull")
plot(x2, pch = 16, cex = 3, border = "red", add = TRUE)

x2 <- st_union(x) |>
    st_concave_hull(ratio = 0.2) |>
    st_geometry()

plot(st_geometry(x), main = "st_concave_hull")
plot(x2, pch = 16, cex = 3, border = "red", add = TRUE)

x2 <- st_buffer(x, dist = 2000) |>
    st_geometry()

plot(st_geometry(x), main = "st_buffer on each point")
plot(x2, pch = 16, cex = 3, border = "red", add = TRUE)

x2 <- st_union(x) |>
    st_buffer(dist = 2000) |>
    st_geometry()

plot(st_geometry(x), main = "st_buffer on all points")
plot(x2, pch = 16, cex = 3, border = "red", add = TRUE)

A note on st_union

The function st_unionis used to combine geometries together. This is useful if we want to apply an operation on each feature or on all combined features.

dim(x)

## [1]  53 131

dim(st_union(x))

## NULL

length(st_union(x))

## [1] 1

head(x)  # several points

## Simple feature collection with 6 features and 130 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -573319.3 ymin: 188315 xmax: -539407.6 ymax: 202856
## Projected CRS: NAD83 / Quebec Lambert
## # A tibble: 6 × 131
##   key        scientificName      issues      datasetKey   publishingOrgKey installationKey
##   <chr>      <chr>               <chr>       <chr>        <chr>            <chr>          
## 1 4102609133 Trillium erectum L. cdc,cdround 50c9509d-22… 28eb1a3f-1c15-4… 997448a8-f762-…
## 2 4103164245 Trillium erectum L. cdc,cdround 50c9509d-22… 28eb1a3f-1c15-4… 997448a8-f762-…
## 3 4103246718 Trillium erectum L. cdc,cdround 50c9509d-22… 28eb1a3f-1c15-4… 997448a8-f762-…
## 4 4103308652 Trillium erectum L. cdc,cdround 50c9509d-22… 28eb1a3f-1c15-4… 997448a8-f762-…
## 5 4111654289 Trillium erectum L. cdc,cdround 50c9509d-22… 28eb1a3f-1c15-4… 997448a8-f762-…
## 6 4112004719 Trillium erectum L. cdc,cdround 50c9509d-22… 28eb1a3f-1c15-4… 997448a8-f762-…
## # ℹ 125 more variables: hostingOrganizationKey <chr>, publishingCountry <chr>,
## #   protocol <chr>, lastCrawled <chr>, lastParsed <chr>, crawlId <int>,
## #   basisOfRecord <chr>, occurrenceStatus <chr>, taxonKey <int>, kingdomKey <int>,
## #   phylumKey <int>, classKey <int>, orderKey <int>, familyKey <int>, genusKey <int>,
## #   speciesKey <int>, acceptedTaxonKey <int>, acceptedScientificName <chr>,
## #   kingdom <chr>, phylum <chr>, order <chr>, family <chr>, genus <chr>, species <chr>,
## #   genericName <chr>, specificEpithet <chr>, taxonRank <chr>, taxonomicStatus <chr>, …

head(st_union(x))  # on set of MULTIPOINT

## Geometry set for 1 feature 
## Geometry type: MULTIPOINT
## Dimension:     XY
## Bounding box:  xmin: -577795.8 ymin: 185325.3 xmax: -538252 ymax: 203526.6
## Projected CRS: NAD83 / Quebec Lambert

Operations on polygons

First, let’s create two polygons with buffers around some Red Trillium observations.

# construct two buffers
x1 <- st_buffer(x[38, ], dist = 4000)
x2 <- st_buffer(x[23, ], dist = 2000)

plot(st_geometry(rbind(x1, x2)), axes = TRUE)  # bind them in the same object

Overlapping polygons

Here are some common overlapping operations that are done on polygons.

# list of functions to apply
op <- c("st_intersection", "st_difference", "st_sym_difference", "st_union")
# set the graphic window
par(mfrow = c(2, 2), mar = c(0, 0, 3, 0))

limits <- st_union(x1, x2) |>
    st_bbox()
xlim <- limits[c(1, 3)]
ylim <- limits[c(2, 4)]

for (i in seq_along(op)) {
    plot(st_geometry(x1), main = op[i], xlim = xlim, ylim = ylim)
    plot(st_geometry(x2), add = TRUE)
    plot(st_geometry(get(op[i])(x1, x2)), add = TRUE, col = "forestgreen")  # this is just a way to get all functions defined in op
}

Measure distances

Here, let’s measure the distance between each points and their centroid.

cx <- st_union(x) |>
    st_centroid()
plot(st_geometry(x), axes = TRUE)
plot(st_geometry(cx), add = TRUE, col = "red", pch = 16, lwd = 3)
text(st_coordinates(x), label = round(as.numeric(st_distance(x, cx))/1000, 1), adj = c(-0.5,
    1.5), cex = 0.6)
mtext(side = 3, line = 0, text = "Distance in km from centroid", ad = 0.98)

Measure areas

set.seed(123)

pol <- st_buffer(x[1:2, ], dist = runif(nrow(x[1:2, ]), 1000, 3000))
area <- st_area(pol)
plot(st_geometry(pol), axes = TRUE)
text(st_coordinates(st_centroid(pol)), label = paste(round(as.numeric(area)/1000), "km"))

Sample

Sample points in polygons in different ways using the st_sample function.

set.seed(12345)
n <- 20
type <- c("random", "regular", "hexagonal")  # list sampling types
par(mfrow = c(2, 2), mar = c(0, 0, 2, 0))  # set graphic window

# plot the different sampling types
for (i in seq_along(type)) {
    s <- st_sample(pol, n, type = type[i])  # get random points
    plot(st_geometry(pol), main = type[i], font = 2)  # plot polygon
    plot(st_geometry(s), pch = 1, add = TRUE)  # plot random points
}
s <- st_sample(pol, rep(n, nrow(pol)), type = type[1])  # get random points
plot(st_geometry(pol), main = type[1], font = 2)  # plot polygon
plot(st_geometry(s), pch = 1, add = TRUE)  # plot random points

Combine operations

Here is an example were we combine several operations to determine the length of roads within each random buffer in a region using the previous road data.

Let’s build the random buffers and see the results.

set.seed(1234)  # set the seed to obtain the same result each time
n <- 10
b <- st_convex_hull(st_union(roads))  # determine the bounding box of the road system
s <- st_sample(b, n, type = "random")  # throw random points within the bounding box
buffs <- st_buffer(s, dist = 10000)  # turn the points to buffers

plot(st_geometry(roads), col = gray(0.7), axes = TRUE)
text(st_coordinates(s), labels = names(buffs))
plot(buffs, border = NA, col = gray(0, 0.1), add = TRUE)

Extract all road sections that intersect with a buffer.

int <- st_intersection(buffs, roads, byid = TRUE)
plot(int, col = "red", lwd = 4, add = TRUE)

Then determine the list of road sections touching each buffer.

o <- st_intersects(buffs, int)
o

## Sparse geometry binary predicate list of length 10, where the predicate was
## `intersects'
##  1: 3, 4, 5, 6
##  2: 14, 15, 16, 17, 45, 70, 76, 77, 94, 95, ...
##  3: 9, 18, 19, 43, 44, 50, 51, 52, 58, 59, ...
##  4: (empty)
##  5: 78, 79, 80, 81, 82, 83, 84, 98, 99, 100, ...
##  6: (empty)
##  7: (empty)
##  8: 7, 8
##  9: 1, 2, 10, 11, 12, 13, 20, 21, 22, 23, ...
##  10: (empty)

Finally, calculate the total length of roads inside each buffer.

dist <- sapply(o, function(i) {
    sum(st_length(int[i]))
})
text(st_coordinates(st_centroid(buffs)), labels = round(dist/1000, 1), col = "black", adj = c(-0.75,
    -1.5), xpd = TRUE, cex = 1)

Simplify things

plot(st_geometry(canada))
plot(st_geometry(st_simplify(canada, preserveTopology = TRUE, dTolerance = 50000)), border = "red",
    add = TRUE)

ms_simplify(canada, 0.2) |>
    st_geometry() |>
    plot()

– Exercice 3 –

Calculate the proportion of forests in 10 random buffers across a landscape. Let’s first generate a fictious landscape.

set.seed(123)

#x<-runif(100, -74, -71) # generate random points
#y<-runif(100, 45, 47)


h <- x |> # turn it to spatial points 
     st_as_sfc() |> 
     st_buffer(dist = 300 * rgamma(nrow(x), 2, 1)) |> # build random gamma buffers around them 
     lapply(function(i){st_as_sfc(st_bbox(i))}) |>
     do.call("c", args = _) |>
     st_union() # put them all in a single polygon
   
plot(h, col = "green4", border = NA, axes = TRUE)

# First throw 10 random points within the envelope of the landscape and then generates buffers around these points
b <- st_convex_hull(h) |>
     st_sample(10, type = "random") |>
     st_buffer(dist = 1000)


# Plot the buffers
plot(h, col = "green4", border = NA, axes = TRUE)
plot(b, add = TRUE)

# Intersect the buffers with the forest and plot the result
cuts <- st_intersection(b, h)
plot(cuts, add = TRUE, col = "red", border = NA)

int <- lengths(st_intersects(b,h))

int[ int > 0] <- st_area(cuts)

perc <- 100 * int / st_area(b)

text(st_coordinates(st_centroid(b)), label = round(perc, 2), cex = 0.7)

RASTERS (terra)

What is a raster?

A raster is a regular grid of pixel with values. Here is an example of building a simple raster with random values using the terra package.

library(terra)

n <- 200
r <- rast(nrow = n, ncol = n, ext = ext(canada), crs = crs(canada))
r <- setValues(r, runif(n^2, 0, 1))
ncell(r)

## [1] 40000

plot(r)
plot(st_geometry(canada), add = TRUE)

Read and write

Writing a raster to a file is done with the writeRaster function.

writeRaster(r, "/home/frousseu/Documents/github/spatialR/data/rast.tif", overwrite = TRUE)
list.files("/home/frousseu/Documents/github/spatialR/data")

##  [1] "carreteras.dbf"        "carreteras.prj"        "carreteras.sbn"       
##  [4] "carreteras.sbx"        "carreteras.shp"        "carreteras.shp.xml"   
##  [7] "carreteras.shx"        "cat.txt"               "cdem_dem_021E_tif"    
## [10] "cdem_dem_021E.tif"     "cdem_dem_031F_tif"     "cdem_dem_031F_tif.zip"
## [13] "cdem_dem_031F.tif"     "cdem_dem_031G_tif"     "cdem_dem_031G_tif.zip"
## [16] "cdem_dem_031G.tif"     "cdem_dem_031H_tif"     "cdem_dem_031H.tif"    
## [19] "rast.tif"              "roads.dbf"             "roads.gpkg"           
## [22] "roads.prj"             "roads.shp"             "roads.shx"            
## [25] "test.dbf"              "test.gpkg"             "test.prj"             
## [28] "test.shp"              "test.shx"

Reading can be done directly with the rast function.

r1 <- rast("/home/frousseu/Documents/github/spatialR/data/cdem_dem_031F.tif")
r2 <- rast("/home/frousseu/Documents/github/spatialR/data/cdem_dem_031G.tif")

par(mfrow = 1:2)
plot(r1)
plot(r2)

Aggregate, merge and crop

Rasters can be aggregated, merged or cropped. An aggregation reduces the number of cells by a certain factor, while a merge brings two rasters together.

ncell(r1)

## [1] 46094401

r1 <- aggregate(r1, fact = 10)
r2 <- aggregate(r2, fact = 10)

ncell(r1)

## [1] 462241

r <- terra::merge(r1, r2)  # mosaic could also be used if rasters overlap
plot(r)

Rasters can be cropped with a specific region using the extent argument.

e <- ext(-76.5, -75, 45.2, 45.8)
plot(crop(r, e))
plot(st_geometry(outaouais), add = TRUE)

plot(crop(r, outaouais))
plot(st_geometry(outaouais), add = TRUE)

crop(r, outaouais) |>
    mask(outaouais) |>
    plot()
plot(st_geometry(outaouais), add = TRUE)

Dimensions

Rasters can have a single or multiple layers. However, layers need to be of the same exact dimensions and resolution.

r

## class       : SpatRaster 
## dimensions  : 481, 1921, 1  (nrow, ncol, nlyr)
## resolution  : 0.002083333, 0.002083333  (x, y)
## extent      : -78.0001, -73.99802, 44.99802, 46.0001  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat NAD83(CSRS) (EPSG:4617) 
## source(s)   : memory
## name        : cdem_dem_031F 
## min value   :          20.0 
## max value   :         588.7

dim(r)

## [1]  481 1921    1

ncell(r)

## [1] 924001

nlyr(r)

## [1] 1

Memory

When reading rasters on your computer, often they will not be brought into memory. The object created will just be a pointer to a file. This is to reduce the memory footprint of rasters.

inMemory(r)

## [1] TRUE

lc <- rast("WorldCover_trees_30s.tif")
inMemory(lc)

## [1] FALSE

Cropping and masking

lc <- rast("WorldCover_trees_30s.tif")
inMemory(lc)

## [1] FALSE

f <- crop(lc, canada)  # crop to Canada extent
plot(f)

f <- mask(f, canada)  # mask to the Canada geometry
plot(f)

f <- crop(lc, canada, mask = TRUE)  # crop and mask to Canada extent
plot(f)

Projecting

f <- project(lc, r)

The function project will take what is in lc and put it in the resolution and extent of r.

Stack rasters together

forest <- f
names(f) <- "forest"

elevation <- r
names(elevation) <- "elevation"

env <- c(forest, elevation)

dim(env)

## [1]  481 1921    2

env

## class       : SpatRaster 
## dimensions  : 481, 1921, 2  (nrow, ncol, nlyr)
## resolution  : 0.002083333, 0.002083333  (x, y)
## extent      : -78.0001, -73.99802, 44.99802, 46.0001  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat NAD83(CSRS) (EPSG:4617) 
## source(s)   : memory
## names       : trees, elevation 
## min values  :     0,      20.0 
## max values  :     1,     588.7

A note on reprojecting rasters

Each time a raster is projected or transformed, the original values will be modified and calculation will be done to fit it in the new dimensions and resolution of the output raster. Original values are lost. Thus, it is always better to transform the crs of vectors to fit it to a raster than vice versa.

– Exercice 4 –

Crop and mask an elevation raster to the Quebec province. You can use the geodata package functions to get the elevation data.

library(geodata)
e <- elevation_30s(country = "CAN", path = getwd(), mask = TRUE)

Raster operations

Extract

Let’s extract elevation data for the Red Trillium occurrences we got previously.

x <- trillein
x <- st_transform(x, st_crs(elevation))

plot(elevation)
plot(st_geometry(outaouais), add = TRUE)
plot(st_geometry(x), add = TRUE)

e <- extract(elevation, x)
head(e)

##   ID elevation
## 1  1    102.70
## 2  2     91.96
## 3  3    119.40
## 4  4    117.97
## 5  5    120.76
## 6  6    139.89

elevation

## class       : SpatRaster 
## dimensions  : 481, 1921, 1  (nrow, ncol, nlyr)
## resolution  : 0.002083333, 0.002083333  (x, y)
## extent      : -78.0001, -73.99802, 44.99802, 46.0001  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat NAD83(CSRS) (EPSG:4617) 
## source(s)   : memory
## name        : elevation 
## min value   :      20.0 
## max value   :     588.7

hist(e[, 2], xlab = "Altitude", main = "")

Extract multiple layers

e <- extract(env, x)
head(e)

##   ID     trees elevation
## 1  1 0.6641760    102.70
## 2  2 0.2906505     91.96
## 3  3 0.9656990    119.40
## 4  4 0.9591519    117.97
## 5  5 0.3617998    120.76
## 6  6 0.9340672    139.89

exactextractr

exactextractr is an R package designed for faster extractions using C++. It is usually faster than the function terra::extract. It also makes exact extractions, i.e. it takes into account the proportion of cells that are overlapping with the polygons. This can also be done with terra::extract and the argument exact = TRUE, but it is much slower.

Reclassify

mn <- c(0, 50, 100, 150, 200, 400)  # min
mx <- c(50, 100, 150, 200, 400, 600)  # max
mat <- cbind(mn, mx, lab = mx)
mat

##       mn  mx lab
## [1,]   0  50  50
## [2,]  50 100 100
## [3,] 100 150 150
## [4,] 150 200 200
## [5,] 200 400 400
## [6,] 400 600 600

elevationclass <- classify(elevation, mat)
plot(elevationclass, col = rev(terrain.colors(nrow(mat))))
plot(st_geometry(outaouais), add = TRUE)

Rasterize

Here is an example of building a raster where each pixel value is determined by the number of points in each cell. For this, we use the rasterize function.

trille <- st_transform(trille, 32198)

rtrille <- rast(resolution = 20000, ext = ext(trille))

rtrille <- rasterize(trille, rtrille, field = 1, fun = "count", background = 0)

rtrille[rtrille == 0] <- NA  # replace cells with 0 observation with NA's to better visualize data

rtrille

## class       : SpatRaster 
## dimensions  : 27, 48, 1  (nrow, ncol, nlyr)
## resolution  : 20000, 20000  (x, y)
## extent      : -812437.4, 147562.6, 121020.9, 661020.9  (xmin, xmax, ymin, ymax)
## coord. ref. :  
## source(s)   : memory
## name        : count 
## min value   :     1 
## max value   :   156

### Plot the results

plot(rtrille)
# plot(st_geometry(trille), col = gray(0, 0.25), add = TRUE)
plot(st_geometry(st_transform(canada, st_crs(trille))), add = TRUE)

Build a rasters from distances

Measure distances from roads and store the values in a raster.

roads2 <- st_transform(roads, 32198)

r <- rast(resolution = 1000, ext = ext(roads2))

rcentroids <- xyFromCell(r, 1:ncell(r)) |>
    as.data.frame() |>
    st_as_sf(coords = c("x", "y"), crs = st_crs(roads2))

g <- st_distance(rcentroids, st_union(roads2))

class(g)

## [1] "units"

dim(g)

## [1] 59024     1

r <- setValues(r, g[, 1])

plot(r)
plot(st_geometry(roads2), add = TRUE)

Rasterize with lengths

An example on using rasterizeGeom to compute the length of roads within each pixel of a raster.

test <- rasterizeGeom(vect(roads2), aggregate(r, 10), fun = "length")
plot(test)
plot(st_geometry(roads2), add = TRUE)

Contours

Convert rasters values to contours in a SpatialLines object.

con <- as.contour(r)
con

##  class       : SpatVector 
##  geometry    : lines 
##  dimensions  : 12, 1  (geometries, attributes)
##  extent      : -2686644, -2439644, -2757371, -2520371  (xmin, xmax, ymin, ymax)
##  coord. ref. :  
##  names       :   level
##  type        :   <num>
##  values      :    5000
##                  1e+04
##                1.5e+04

plot(r)
plot(con, add = TRUE)

Polygonize a raster

Turn all cells that are at least 25km from the roads to polygons.

r25 <- r
r25[r25 < 25000] <- NA

pr <- as.polygons(r25, aggregate = TRUE) |>
    st_as_sf()
plot(r)
plot(st_geometry(pr), col = "blue", border = NA, add = TRUE)

– Exercice 5 –

Get data for Trillium cernuum using package rgbif and calculate the mean annual temperature in the species range. Bioclimatic variables are described in the WorldClim pages.

Here is how to download the data. Note that the object returned by occ_search is a gbif object from which the data element need to be extracted. This element is a tibble that needs to be turned to a data.frame in order to convert it to a spatial object.

library(geodata)
library(rgbif)

# download temperature data
tavg <- worldclim_country("CAN", var = "tavg", path = getwd())
tavg <- mean(tavg)

# download records
x <- occ_search(scientificName = "Trillium erectum", limit = 1000, hasCoordinate = TRUE, country = "CA")
x <- as.data.frame(x$data)

Mapping

ggplot2

library(ggplot2)
library(tidyterra)


ggplot() + geom_spatraster(data = tavg) + geom_sf(data = trille, color = "orange", alpha = 0.8) +
    geom_sf(data = canada, color = "white", fill = NA, alpha = 0.8) + coord_sf(xlim = ext(tavg)[1:2],
    ylim = ext(tavg)[3:4], expand = FALSE) + theme_light()

leaflet

library(leaflet)
library(leaflet.extras)

trille <- st_transform(trille, 4326)

leaflet(cat) %>% 
  addTiles() %>%
  addCircleMarkers(data = trille, group = 'trille') %>%
  addDrawToolbar(targetGroup = 'trille', editOptions = editToolbarOptions(selectedPathOptions = selectedPathOptions())) %>%
  addLayersControl(overlayGroups = c('trille'),options =layersControlOptions(collapsed=FALSE)) %>%
  addMeasurePathToolbar(options =measurePathOptions(imperial = TRUE,minPixelDistance = 100,showDistances = FALSE))

mapedit

Create shapes with mapedit functionality.

library(mapview)
library(mapedit)

x <- mapview(trille) |>
    editMap()

x <- x$drawn

plot(st_geometry(x))

mapview(trille) + mapview(x)

Other resources

rspatial.org

Geocomputation with R